# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating figures and tables in Appendix G
# Appendix G Table 21: Influence of Racial Boundary on Pedestrian Stops
# Appendix G Figure 12: Influence of Racial Boundary on Police Stops (Standardized). Annotations denote the robustness value and bounding variable value necessary to attenuate the substantive influence of racial boundaries to 0.
# Appendix G Table 22: Influence of Logged Crime on Logged Pedestrian Stops (Standardized), Conditional on Racial Boundary Status.
# Appendix G Figure 13: Influence of Logged Crime on Pedestrian Stops, Conditional on Racial Boundary Status.
# Appendix G Table 23: Influence of Percent White on Logged Pedestrian Stops (Standardized), Conditional on Racial Boundary Status
# Appendix G Figure 14: Influence of Percent White on Pedestrian Stops, Conditional on Racial Boundary Status.
# Table 24: Influence of Percent White on Pedestrian Stops by Race of Civilian, Conditional on Racial Boundary Status
# Figure 15: Influence of Percent White on Pedestrian Stops by Race of Civilian, Conditional on Racial Boundary Status.


suppressPackageStartupMessages(
  
  {
    library(AER)
    library(dplyr)
    library(fixest)
    library(lfe)
    library(tidyverse)
    library(ggplot2)
    library(MASS)
    library(sensemakr)
    library(haven)
    library(readstata13)
    library(readxl)
    library(readr)
    library(gridExtra)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)


# load stop data
load("aus_stops_final.RData")
# now run analyses

# log and +1 to variables (pop already logged)
aus_stops_final = aus_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
aus_stops_final = aus_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1),
         lped_stops_total = log(ped_stops_total + 1),
         lped_stops_black = log(ped_stops_black + 1),
         lped_stops_latino = log(ped_stops_latino + 1),
         lped_stops_white = log(ped_stops_white + 1),
         lped_stops_asian = log(ped_stops_asian + 1),
         lped_stops_nonwhite = log(ped_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_stops_final$larrest_sd <- aus_stops_final$larrests - (mean(aus_stops_final$larrests)/sd(aus_stops_final$larrests))
aus_stops_final$lmisdemeanors_sd <- aus_stops_final$lmisdemeanors - (mean(aus_stops_final$lmisdemeanors)/sd(aus_stops_final$lmisdemeanors))
aus_stops_final$lfelonies_sd <- aus_stops_final$lfelonies - (mean(aus_stops_final$lfelonies)/sd(aus_stops_final$lfelonies))
aus_stops_final$lnonviolent_sd <- aus_stops_final$lnonviolent - (mean(aus_stops_final$lnonviolent)/sd(aus_stops_final$lnonviolent))
aus_stops_final$lviolent_sd <- aus_stops_final$lviolent - (mean(aus_stops_final$lviolent)/sd(aus_stops_final$lviolent))
aus_stops_final$lsociety_sd <- aus_stops_final$lsociety - (mean(aus_stops_final$lsociety)/sd(aus_stops_final$lsociety))
aus_stops_final$lperson_sd <- aus_stops_final$lperson - (mean(aus_stops_final$lperson)/sd(aus_stops_final$lperson))
aus_stops_final$lproperty_sd <- aus_stops_final$lproperty - (mean(aus_stops_final$lproperty)/sd(aus_stops_final$lproperty))

# create scaled DVs for stop vars
aus_stops_final$lall_stops_total_sd <- aus_stops_final$lall_stops_total - (mean(aus_stops_final$lall_stops_total)/sd(aus_stops_final$lall_stops_total))
aus_stops_final$lall_stops_black_sd <- aus_stops_final$lall_stops_black - (mean(aus_stops_final$lall_stops_black)/sd(aus_stops_final$lall_stops_black))
aus_stops_final$lall_stops_latino_sd <- aus_stops_final$lall_stops_latino - (mean(aus_stops_final$lall_stops_latino)/sd(aus_stops_final$lall_stops_latino))
aus_stops_final$lall_stops_white_sd <- aus_stops_final$lall_stops_white - (mean(aus_stops_final$lall_stops_white)/sd(aus_stops_final$lall_stops_white))
aus_stops_final$lall_stops_asian_sd <- aus_stops_final$lall_stops_asian - (mean(aus_stops_final$lall_stops_asian)/sd(aus_stops_final$lall_stops_asian))
aus_stops_final$lall_stops_nonwhite_sd <- aus_stops_final$lall_stops_nonwhite - (mean(aus_stops_final$lall_stops_nonwhite)/sd(aus_stops_final$lall_stops_nonwhite))

# create binned boundary measure
aus_stops_final <- aus_stops_final %>% mutate(boundary_quart = quantile(aus_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# Milwaukee
load("mil_stops_final.RData")


# now run analyses

# log and +1 to variables (pop already logged)
mil_stops_final = mil_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
mil_stops_final = mil_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_stops_final$larrest_sd <- mil_stops_final$larrests - (mean(mil_stops_final$larrests)/sd(mil_stops_final$larrests))
mil_stops_final$lmisdemeanors_sd <- mil_stops_final$lmisdemeanors - (mean(mil_stops_final$lmisdemeanors)/sd(mil_stops_final$lmisdemeanors))
mil_stops_final$lfelonies_sd <- mil_stops_final$lfelonies - (mean(mil_stops_final$lfelonies)/sd(mil_stops_final$lfelonies))
mil_stops_final$lnonviolent_sd <- mil_stops_final$lnonviolent - (mean(mil_stops_final$lnonviolent)/sd(mil_stops_final$lnonviolent))
mil_stops_final$lviolent_sd <- mil_stops_final$lviolent - (mean(mil_stops_final$lviolent)/sd(mil_stops_final$lviolent))
mil_stops_final$lsociety_sd <- mil_stops_final$lsociety - (mean(mil_stops_final$lsociety)/sd(mil_stops_final$lsociety))
mil_stops_final$lperson_sd <- mil_stops_final$lperson - (mean(mil_stops_final$lperson)/sd(mil_stops_final$lperson))
mil_stops_final$lproperty_sd <- mil_stops_final$lproperty - (mean(mil_stops_final$lproperty)/sd(mil_stops_final$lproperty))

# create scaled DVs for stop vars
mil_stops_final$lall_stops_total_sd <- mil_stops_final$lall_stops_total - (mean(mil_stops_final$lall_stops_total)/sd(mil_stops_final$lall_stops_total))
mil_stops_final$lall_stops_black_sd <- mil_stops_final$lall_stops_black - (mean(mil_stops_final$lall_stops_black)/sd(mil_stops_final$lall_stops_black))
mil_stops_final$lall_stops_latino_sd <- mil_stops_final$lall_stops_latino - (mean(mil_stops_final$lall_stops_latino)/sd(mil_stops_final$lall_stops_latino))
mil_stops_final$lall_stops_white_sd <- mil_stops_final$lall_stops_white - (mean(mil_stops_final$lall_stops_white)/sd(mil_stops_final$lall_stops_white))
mil_stops_final$lall_stops_asian_sd <- mil_stops_final$lall_stops_asian - (mean(mil_stops_final$lall_stops_asian)/sd(mil_stops_final$lall_stops_asian))
mil_stops_final$lall_stops_nonwhite_sd <- mil_stops_final$lall_stops_nonwhite - (mean(mil_stops_final$lall_stops_nonwhite)/sd(mil_stops_final$lall_stops_nonwhite))

# create binned boundary measure
mil_stops_final <- mil_stops_final %>% mutate(boundary_quart = quantile(mil_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

# Chicago
load("chi_stops_final.RData")

# now run analyses

# log and +1 to variables (pop already logged)
chi_stops_final = chi_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime + 1),
         lviolentcrime = log(violent_crime + 1))

# now for new stop vars
chi_stops_final = chi_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_stops_final$larrest_sd <- chi_stops_final$larrests - (mean(chi_stops_final$larrests)/sd(chi_stops_final$larrests))
chi_stops_final$lmisdemeanors_sd <- chi_stops_final$lmisdemeanors - (mean(chi_stops_final$lmisdemeanors)/sd(chi_stops_final$lmisdemeanors))
chi_stops_final$lfelonies_sd <- chi_stops_final$lfelonies - (mean(chi_stops_final$lfelonies)/sd(chi_stops_final$lfelonies))
chi_stops_final$lnonviolent_sd <- chi_stops_final$lnonviolent - (mean(chi_stops_final$lnonviolent)/sd(chi_stops_final$lnonviolent))
chi_stops_final$lviolent_sd <- chi_stops_final$lviolent - (mean(chi_stops_final$lviolent)/sd(chi_stops_final$lviolent))
chi_stops_final$lsociety_sd <- chi_stops_final$lsociety - (mean(chi_stops_final$lsociety)/sd(chi_stops_final$lsociety))
chi_stops_final$lperson_sd <- chi_stops_final$lperson - (mean(chi_stops_final$lperson)/sd(chi_stops_final$lperson))
chi_stops_final$lproperty_sd <- chi_stops_final$lproperty - (mean(chi_stops_final$lproperty)/sd(chi_stops_final$lproperty))

# create scaled DVs for stop vars
chi_stops_final$lall_stops_total_sd <- chi_stops_final$lall_stops_total - (mean(chi_stops_final$lall_stops_total)/sd(chi_stops_final$lall_stops_total))
chi_stops_final$lall_stops_black_sd <- chi_stops_final$lall_stops_black - (mean(chi_stops_final$lall_stops_black)/sd(chi_stops_final$lall_stops_black))
chi_stops_final$lall_stops_latino_sd <- chi_stops_final$lall_stops_latino - (mean(chi_stops_final$lall_stops_latino)/sd(chi_stops_final$lall_stops_latino))
chi_stops_final$lall_stops_white_sd <- chi_stops_final$lall_stops_white - (mean(chi_stops_final$lall_stops_white)/sd(chi_stops_final$lall_stops_white))
chi_stops_final$lall_stops_asian_sd <- chi_stops_final$lall_stops_asian - (mean(chi_stops_final$lall_stops_asian)/sd(chi_stops_final$lall_stops_asian))
chi_stops_final$lall_stops_nonwhite_sd <- chi_stops_final$lall_stops_nonwhite - (mean(chi_stops_final$lall_stops_nonwhite)/sd(chi_stops_final$lall_stops_nonwhite))

# create binned boundary measure
chi_stops_final <- chi_stops_final %>% mutate(boundary_quart = quantile(chi_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# Austin, DV = total stops

# rename racial and economic boundary variable
aus_stops_final$aus_white_blv <- aus_stops_final$p_race_white_blv
aus_stops_final$aus_ses_blv <- aus_stops_final$ses_blv



#### reg tables associated with boundary --> ped stops ####

aus.blv.ped.stops = lm_robust(lped_stops_total_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, 
                              data = aus_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.blv.ped.stops = lm_robust(lped_stops_total_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, 
                              data = chi_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for white_blv --> stops
texreg(l = list(aus.blv.ped.stops, chi.blv.ped.stops),
       file = "regtable.ped.stopsmain.tex",
       include.ci = FALSE,
       custom.coef.map = list("p_race_white_blv" = "Boundary (White)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Racial Boundary on Pedestrian Stops",
       custom.header = list("Dependent Variable: Logged Pedestrian Stops" = 1:2),
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("AUS","CHI"),
       label = "table:ped_stops_mainfx",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


## white boundary measure on all stops
aus.model.ped.stops = lm_robust(lped_stops_total_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                                data = aus_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(aus.model.ped.stops)

aus.model.ped.stops.sense = lm(lped_stops_total_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                               data = aus_stops_final, subset = pop > 0)


## create df of austin coefficients
coef.aus.1 <- tidy(aus.model.ped.stops)

fit_cis_95 <- confint(aus.model.ped.stops, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(aus.model.ped.stops, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.aus.1 <- bind_cols(coef.aus.1, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
aus_boundary_white <- coef.aus.1[1,]


# Chicago DV = ped stops
# rename racial and economic boundary variable
chi_stops_final$chi_white_blv <- chi_stops_final$p_race_white_blv
chi_stops_final$chi_ses_blv <- chi_stops_final$ses_blv

## white boundary measure on all stops
chi.model.ped.stops = lm_robust(lped_stops_total_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                                data = chi_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(chi.model.ped.stops)

chi.model.ped.stops.sense = lm(lped_stops_total_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                               data = chi_stops_final, subset = pop > 0)


## create df of chicago coefficients
coef.chi.1 <- tidy(chi.model.ped.stops)

fit_cis_95 <- confint(chi.model.ped.stops, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(chi.model.ped.stops, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.chi.1 <- bind_cols(coef.chi.1, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
chi_boundary_white <- coef.chi.1[1,]


## combine 2 dfs w white boundary coefs
combined.coeffs <- rbind(aus_boundary_white, chi_boundary_white)

combined.coeffs <- combined.coeffs %>% dplyr::select(-c(SE,statistic,p.value, conf.low, conf.high))

## make xaxis labels
Plot1labels <- c("Austin","Chicago")

## sensitivity analysis

# robustness value exercise 
library(sensemakr)
sense_out1 = 
  sensemakr(aus.model.ped.stops.sense,
            treatment = "aus_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

sense_out2 = 
  sensemakr(chi.model.ped.stops.sense,
            treatment = "chi_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 35)


rvdf = data.frame(
  
  rv_val = c(as.numeric(sense_out1$sensitivity_stats$rv_q),
             as.numeric(sense_out2$sensitivity_stats$rv_q)),
  
  variable = combined.coeffs$Variable,
  
  bound = c("8x Violent Crime", "35x Violent Crime")
  
)


combined.coeffs$se = (combined.coeffs$conf.high_95 - combined.coeffs$Coefficient) / 1.96
# combined.coeffs = combined.coeffs %>% filter(Variable != 'balt_white_blv')

library(meta)
out_meta = metagen(
  TE = combined.coeffs$Coefficient,
  seTE = combined.coeffs$se
)

# Austin: 9x violent crime
# Chicago: 35x violent crime


df_meta = data.frame(
  Variable = "Meta-analysis",
  Coefficient = out_meta$TE.random,
  df = out_meta$df.Q,
  outcome = "lped_stops_total_sd",
  se = out_meta$seTE.random,
  conf.low_95 = out_meta$TE.random - 1.96 * out_meta$seTE.random,
  conf.high_95 = out_meta$TE.random + 1.96 * out_meta$seTE.random,
  conf.low_90 = out_meta$TE.random - 1.645 * out_meta$seTE.random,
  conf.high_90 = out_meta$TE.random + 1.645 * out_meta$seTE.random
)

combined.coeffs = bind_rows(combined.coeffs, df_meta)
combined.coeffs$Variable = 
  factor(combined.coeffs$Variable, levels = combined.coeffs$Variable)

library(ggthemes)
white_blv_plot_pedstops <- 
  ggplot(combined.coeffs, 
         aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 1) +
  geom_point(aes(x = Variable, 
                 y = Coefficient)) + 
  geom_errorbar(aes(x = Variable,
                    ymin = conf.low_90,
                    ymax = conf.high_90),
                width = 0,
                size = .8) +
  geom_errorbar(aes(x = Variable,
                    ymin = conf.low_95,
                    ymax = conf.high_95),
                width = 0,
                size = .4) +
  annotate("text",
           label =  paste0("RV: ", round(rvdf$rv_val[1:2], 2),
                           "\n", rvdf$bound[1:2]),
           y = combined.coeffs$Coefficient[1:2],
           x = c(1.4, 2.4),
           size = 4,
           family = "Times") +
  geom_hline(yintercept = combined.coeffs$Coefficient[3], 
             colour = gray(1/2), lty = 2,
             alpha = .4,
             size = .2) + 
  xlab("Cities") + ylab("Racial Boundary Coefficients") + 
  scale_color_grey(start = .6, end = 0) + 
  theme_tufte(base_size = 9) + 
  scale_x_discrete(labels = c(Plot1labels, "Meta-Analysis")) + 
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(
    plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 18)
  )

ggsave(width = 8, height = 6, plot = white_blv_plot_pedstops, filename = "all_cities_coef_plot_ped_stops.png")


#### reg tables associated with boundaryXcrime --> ped stops ####

aus.ped.stops.crime = lm_robust(lped_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = aus_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.ped.stops.crime = lm_robust(lped_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                                data = chi_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundaryXcrime --> ped stops
texreg(l = list(aus.ped.stops.crime, chi.ped.stops.crime),
       file = "regtable.ped.stops.crime.mod.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:lcrime" = "Boundary X Crime",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Logged Crime on Logged Pedestrian Stops (Standardized), Conditional on Racial Boundary Status.",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("AUS","CHI"),
       custom.header = list("Dependent Variable: Logged Pedestrian Stops" = 1:2),
       label = "table:ped_stops_crimeXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


###### race_blvXcrime plots #####
## Austin

library(prediction)

mdata <- aus_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy,
                                          lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                          hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                          phys_bound, com_density) %>% na.omit()


aus1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(aus1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

aus_stops_preds <- preds %>% mutate(city = "Austin")


## Milwaukee

library(prediction)


mdata <- mil_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy,
                                          lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                          hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                          phys_bound, com_density) %>% na.omit()


mil1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(mil1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

mil_stops_preds <- preds %>% mutate(city = "Milwaukee")


## Chicago

library(prediction)


mdata <- chi_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy,
                                          lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                          hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                          phys_bound, com_density) %>% na.omit()


chi1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,9.3,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(chi1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,9.3,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

chi_stops_preds <- preds %>% mutate(city = "Chicago")


# graph out of sample predictions
aus_stops_preds <- aus_stops_preds %>% dplyr::select(lcrime,fitted.fit,
                                                     fitted.lwr,fitted.upr,boundary,city)

mil_stops_preds <- mil_stops_preds %>% dplyr::select(lcrime,fitted.fit,
                                                     fitted.lwr,fitted.upr,boundary,city)

chi_stops_preds <- chi_stops_preds %>% dplyr::select(lcrime,fitted.fit,
                                                     fitted.lwr,fitted.upr,boundary,city)



d <- rbind(aus_stops_preds, chi_stops_preds, mil_stops_preds)


# make plot
p <- ggplot(data=d, aes(x=lcrime, y=fitted.fit,group=boundary)) +
  facet_wrap(~ city, nrow = 3, scales = "free", strip.position = "right")+
  #geom_point(aes(shape = Boundary), color = "grey",alpha=.5)+
  geom_line(aes(y=fitted.fit, linetype=boundary, color = boundary))+
  geom_ribbon(aes(ymin=fitted.lwr, ymax=fitted.upr,linetype=boundary, color=boundary),alpha=0.1)+
  labs(y="Pedestrian Stops (Logged)",x="Crime (Logged)") +
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(strip.background = element_rect(fill="white"))

ggsave(plot = p, filename = "stops_crime_mod_plots.png", 
       width = 8, height = 12, bg = 'white')


#### reg tables associated with boundary X %white --> ped stops ####

aus.ped.stops.race.mod = lm_robust(lped_stops_total_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                     pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                     phys_bound + com_density, 
                                   data = aus_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.ped.stops.race.mod = lm_robust(lped_stops_total_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
                                     pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                     phys_bound + com_density, 
                                   data = chi_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# Not enough ped stops to run analysis with milwaukee
# mil.stops.race.mod = lm_robust(lped_stops_total_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
#                                  pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
#                                data = mil_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


# create latex table for boundary X %white --> stops
texreg(l = list(aus.ped.stops.race.mod, chi.ped.stops.race.mod),
       file = "regtable.ped.stops.race.mod.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:p_race_white" = "Boundary X % White",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Percent White on Logged Pedestrian Stops (Standardized), Conditional on Racial Boundary Status",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("AUS","CHI"),
       custom.header = list("Dependent Variable: Logged Pedestrian Stops" = 1:2),
       label = "table:ped_stops_whiteXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


#################################################
# Create Plot for boundaryX%white for ped stops #
#################################################

# plot interaction models for all stops - now using binned racial boundary measure
form_list = 
  list(
    'lped_stops_total_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density' )

hetdf_list = list(form_list)
hetdf_lab = c("Pedestrian Stops")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), aus_stops_final)
  
  fakedata = 
    aus_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Austin")


# create plot
aus_ped_stops_race_mod_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x=element_blank())

# ggsave(width = 8, height = 6, plot = aus_allstops_racial_int, filename = "austin/aus_allstops_racial_int_binned.png")


## Chicago

# plot interaction models for ped stops - now using binned racial boundary measure
form_list = 
  list(
    'lped_stops_total_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density' )

hetdf_list = list(form_list)
hetdf_lab = c("Pedestrian Stops")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), chi_stops_final)
  
  fakedata = 
    chi_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Chicago")


# create plot
chi_ped_stops_race_mod_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots



# put all plots together
library(cowplot)
ped_stops_race_mod_plots <- plot_grid(aus_ped_stops_race_mod_plot, chi_ped_stops_race_mod_plot,
                                      align = "v", nrow = 2, rel_heights = c(1/5, 1/5))



ggsave(plot = ped_stops_race_mod_plots, filename = "ped_stops_race_mod_plots.png", 
       width = 8, height = 12, bg = 'white')


## FOW for ped stops ##
#### reg tables associated with boundary X %white --> ped stops by race ####

# Austin
aus.ped.stops.nonwhite = lm_robust(lped_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + 
                                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                     lviolentcrime + phys_bound + com_density, 
                                   data = aus_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus.ped.stops.white = lm_robust(lped_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                  lviolentcrime + phys_bound + com_density, 
                                data = aus_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# Chicago
chi.ped.stops.nonwhite = lm_robust(lped_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + 
                                     hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                     lviolentcrime + phys_bound + com_density, 
                                   data = chi_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.ped.stops.white = lm_robust(lped_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + 
                                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                  lviolentcrime + phys_bound + com_density, 
                                data = chi_stops_final, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundary X %white --> stops by race
texreg(l = list(aus.ped.stops.nonwhite, aus.ped.stops.white, chi.ped.stops.nonwhite, chi.ped.stops.white),
       file = "regtable.ped.stops.fow.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:p_race_white" = "Boundary * % White",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Percent White on Pedestrian Stops by Race of Civilian, Conditional on Racial Boundary Status",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("Non-white","White","Non-white","White"),
       custom.header = list("Dependent Variable: Logged Pedestrian Stops" = 1:4),
       label = "table:ped_stops_FOW",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))



###### race_blvXp_white plots w race of person stopped #####
## Austin

form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Non-white Persons","White Persons")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), aus_stops_final)
  
  fakedata = 
    aus_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Austin")

# create plot
aus_stops_fow_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x=element_blank())


## Chicago
form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Non-white Persons","White Persons")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), chi_stops_final)
  
  fakedata = 
    chi_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Chicago")

# create plot
chi_stops_fow_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots




## Milwaukee
form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Non-white Persons","White Persons")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), mil_stops_final)
  
  fakedata = 
    mil_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Milwaukee")

# create plot
mil_stops_fow_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.text.x = element_blank())


# put all the plots together
library(cowplot)
stops_fow_plots <- plot_grid(aus_stops_fow_plot, chi_stops_fow_plot, mil_stops_fow_plot,
                             align = "v", nrow = 3, rel_heights = c(1/3, 1/3, 1/2.5))

ggsave(plot = stops_fow_plots, filename = "stops_fow_plots.png", 
       width = 8, height = 12, bg = 'white')



